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Abstract. We provide a new, analytic kernel for scale space filtering of 
dMRI data. The kernel is an approximation for the Green’s function of 
a hypo-elliptic diffusion on the 3D rigid body motion group SE(3), for 
fiber enhancement in dMRI. The enhancements are described by linear 
scale space PDEs in the coupled space of positions and orientations em¬ 
bedded in SE(3). As initial condition for the evolution we use either a 
Fiber Orientation Distribution (FOD) or an Orientation Density Func¬ 
tion (ODE). Explicit formulas for the exact kernel do not exist. Although 
approximations well-suited for fast implementation have been proposed 
in literature, they lack important symmetries of the exact kernel. We in¬ 
troduce techniques to include these symmetries in approximations based 
on the logarithm on SE(3), resulting in an improved kernel. Regard¬ 
ing neuroimaging applications, we apply our enhancement kernel (a) to 
improve dMRI tractography results and (b) to quantify coherence of ob¬ 
tained streamline bundles. 

Keywords: Scale Space on SE(3), Contextual enhancement. Left-invariant 
Diffusion, Group convolution, Tractography. 


1 Introduction 

In dMRI it is assumed that axons grouped in bundles in human brain tissue 
restrict the Brownian motion of water molecules such that more diffusion occurs 
along the bundles [T] . By measuring the decay of signal due to diffusion in many 
directions it is possible to obtain information about the underlying microstruc¬ 
ture of the brain tissue and further processing of this data provides clues about 
the anatomical brain connectivity. After a pre-processing procedure in which the 
raw data is corrected for e.g. distortions and motion artefacts, different models 
can be used for further processing of the data [5]. We construct from the data 
a fiber orientation distribution (FOD) function on positions and orientations, 

* The research leading to the results of this article has received funding from the 
European Research Council under the ECs 7th Framework Programme (FP7/2007 
2014)/ERG grant agreement No. 335555. 
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Fig. 1. The pipeline used in this paper for processing dMRI data. CSD provides an 
FOD, that can be enhanced by convolution with our proposed kernel as in Section [3^. 
We visualize the FOD D by a field of glyphs, where at each position y on a grid and 
each orientation n the radius is of the glyph is directly proportional to D(y,n). The 
color corresponds to direction of the orientation. Then, deterministic or probabilistic 
tractography produces fiber bundles, that often contain deviating fibers. They can be 
classified with our proposed coherence quantification, see Section [HJo. 


representing the probability density of finding a fiber in a certain position and 
orientation. For this we use Constrained Spherical Deconvolution (CSD) [3], but 
the methods in this paper can be combined with any model that outputs an 
FOD or an Orientation Density Function (ODF) of water molecules. 

For regularization of dMRI data, various methods exist that include contex¬ 
tual information on position and/or orientation space |4l5l6l7l8l9ll0lllll2j . In 
this paper we pursue the method of scale spaces on the group of positions and 
rotations SE{3). For this we consider diffusions described by the Fokker-Planck 
equations of hypo-elliptic Brownian motion [51 Sect. 4.2]. The effect of these evo¬ 
lutions on the FOD is that elongated structures are enhanced while crossings of 
bundles are preserved. No explicit formulas exist for the exact Green’s function 
of the PDF, but approximations exist in literature [5]: one is the product of two 
SE(2)-kernels and one is based on the logarithmic modulus on SE{3). However, 
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two important invariance properties of the exact kernel are not automatically 
obeyed in these approximations, as pointed out in Section 12.31 

The first (and main) contribution of this paper is that we provide a new an¬ 
alytic kernel approximation for Brownian motion on SE{3) that is, in contrast 
to previous analytic approximations, well-defined on the quotient. As such, it 
carries the appropriate symmetry. Furthermore, the novel, more precise approx¬ 
imation is still well-suited for fast kernel implementations m- 

The second contribution in this work is the application of the improved ap¬ 
proximations in two different places (a) and (b) in the dMRI pipeline as in Fig.[TJ 
We provide two clinically relevant experiments to illustrate this. We apply the 
enhancements to the FOD and show the advantages in the tractography result 
(a) . Secondly, we use the enhancement kernel to construct a density on the space 
of positions and orientations, based on a set of fibers (b). From this we derive a 
measure for the coherence of fibers within a fiber bundle. The symmetry included 
in the new kernel leads to a reduction in computation time of the measure. 

Section [5] covers the theory of the paper and compares previous and new 
approximation kernels. In Section [3] the applications (a) and (b) are presented. 

2 Approximations of Scale Spaces on SE(3) 

First we explain in Section 12.II how the space of 3D positions and orientations, on 
which FODs are defined, can be embedded in the group of 3D rigid body motions 
SE{3). This is followed by a brief introduction to scale spaces on SE{3), see 
Section 12.21 together with a discussion of the two required invariances. In Section 
I2.3l we give two known kernel approximations and we propose an adaptation for 
the logarithmic approximation, such that the desired invariances are induced. 

2.1 The Embedding of x into SE{3) 

Contextual enhancement of a function U : x 5^ —>• R"*", representing an FOD, 

means improving the alignment of elongated structures present in the FOD. 
Such alignment can only be done in a space where positions and orientations are 
coupled. Therefore we embed R^ x in the group of 3D rigid body motions 
SE{3), with group product (x, R)(x',R') = (x -|- Rx', RR'), where R, R' G 
SO{3), the 3D rotation group, and x, x' G R^. Square integrable functions U 
relate to a specific set of functions C/ : SE{3) —>■ R+ via 

U ^ { D(y,n) = D(y,R„) and C/(y,R) = C/(y,Re^) } , (I) 

where from now on, Rn denotes any rotation such that Rne^ = n, and where 
= (0, 0,1)^ denotes our reference axis. The functions U carry a symmetry, 
as they are right-invariant w.r.t. right-action R i— RRe^^a, where from now 
on Ra,i/) denotes the counter-clockwise rotation around axis a by angle ip. This 
implies we should consider left cosets on SE(3), with equivalence relation 

9 ■= (y, R) = g' ■= (y', R') 9~^9' G R y = y' and R' = RRe^.a, (2) 
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for some a G [0, 27r), where from now on H denotes the subgroup of SE(3) whose 
elements are equal to h = The total set of functions U, mentioned 

earlier, is given by: 

h^{SE{3)) := {U G ]L2{SE{3)) \ U{gh) = U{g) for all g G SE{3),hG H}. 

From now on, we identify the function [/ : x —>■ R with a function on the 

group quotient R^ x S'^ := SE{3)/H = S'i?(3)/({0} x SO{2)). Throughout the 
paper we use a tilde to distinguish functions on the group from functions on the 
quotient. Next, we present scale spaces for such functions U and U. 


2.2 Scale Spaces and Group Convolution on R® x and SE(3) 

We define the enhancement evolutions on SE{3) in terms of the left-invariant 
vector fields. They can be considered as differential operators on locally defined 
smooth functions M- Left-invariant vector fields are obtained from a Lie-algebra 
basis for the tangent space Te{SE{3)) at unity element e = (0,/) G SE{3), say 
{Ai ,..., Aq}. This is done with the pushforward (Lg)^ of the left-multiplication 
LgQ = gq- AgU = Ae{U o Lg) = {{Lg)t,Ae){U), for all smooth tj : SE{3) —>■ R. 
Explicit formulas for the vector fields g i—> Ai\g can be obtained by 


AAnU = lim 

t-5-O 


U{gA^^)-U(g) 


j = 1,...,6, 


(3) 


where Te{SE{3)) 9 A i-^ G SE{3) denotes the exponential map yielding the 
1-parameter subgroup | t G R}. Note that Ai\e = Ai. We write 
where {Ali,.^ 2 ,. 4 . 3 } are spatial and {Ai,A 5 ,Ae} are rotation vector fields, for 
which we use the explicit expressions in two different Euler angle parametriza- 
tions as given in [^. Eor geometric intuition, see Fig. [5J Now the (hypo-elliptic) 
diffusion process describing Brownian motion on SE{3) [B], for U on SE{3) can 
be formulated as: 


^W(y,R,t) = D33(Al3)^lT(y,R,t) + D^iAl + Al)W(y,R,t), 

W(y,R,0) = C/(y,R). 


This yields a scale space representation W of the function U with scale parameter 
t > 0. Parameters ZI 33 > 0 and II 44 > 0 influence the amount of spatial and 
angular diffusion, respectively. Solutions for W can be found via finite difference 
approximations of the PDE or via convolution U with a kernel. Fast kernel 
implementations exist m and are particularly suitable for our applications. The 
S'E(3)-convolution with a probability kernel pt ■ SE{3) —?> R+ is given by 

W{g,t) = {Pt*sEi3)U){g) = [ pt{q~^g)U{q)dg. (5) 

JSEi3) 

where dq = dxdR, representing the Haar measure, for all q = (x, i?) G SE{3). 
This evolution is well-defined also on R^ x 5^ and can be written as 

/ = T' 33 (n • V) 2 W(y,n,t) + D44ALBW{y,n,t), 

\W(y,n,0) = C/(y,n), 


(6) 
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Fig. 2. Left: left-invariant vector fields {Ai}i^i in 515(3), restricted to (y,Rn)- Right: 
glyph visualization of the impulse response pt as the scale space kernel of with 
-D 33 = 1 , FI 44 = f?s 5 = 0.02, t = 2. 


with Alb the Laplace-Beltrami operator on 5'^. Solutions of ([5]) are found by 
XI S'^-convolution with the exact solution kernel pt : xi 5^ —>■ R+: 

W^(y,n,t) = (pt U) (y,n) 

= / /pt(R?:,(y-y'),R^-n)-C/(y',n')dy'da(n'). (7) 

52 R 3 

The kernels pt and pt should satisfy certain symmetries, as shown in the following 
definition, lemma and corollary. 

Definition 1. An operator 3 : 'L§{SE{i)) —>• L^(5i5(3)) is called legal if 

{CgO^ = ^oCg 'ig&SE{Z), 

\^oRh = ^ = Rh'oi> yh,h' eH, 

with Cg and Rh the left- and right-regular action on SE(3) respectively, see 

Legal operators 3 induce well-posed operators ^ on L 2 (K^ xi5^) via <P{U){y, n) = 
<^(C/)(y, Rn). In particular, this holds for legal scale space operators. 

Lemma 1. Let ^ : L^(5i5(3)) —>■ L^(5i?(3)) be linear and legal, and assume it 
maps h2{SE{3)) into Loo(S'i?(3)). Then we have: 
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1. identity: {3{U)){g) = J K{g,q)U{q) dq = {p *se{ 3 ) U){g), with 

SE{3) 

Kig, q) = k{e, {g~^q)) =: P{q~^g) and with p{g) := k{g, e) = p{g~^). 

2. symmetry: K{gh, qh') = K{g, q) for all h^h' € H and all g G SE{3). 

3. preservation of correspondences: {U U) ^{U) ^^(U) with 


i^{U)){y,n) = J k(y,n,y',n')U(y',n')dy'da(n'), with 
R3xS2 

k{y, n, y', n') = ^k{{y, R^), (y', Rn')) 

= ^p((y',R-n)”^(y:R-n)) =p(R^/(y-y'),R-r:/n). 

Proof. By the Dunford-Pettis Theorem cf. [TB], is a kernel operator and ||^|P = 
sup / 1 ^( 5 , (?)|2dg, so we can write (^(Lf))(g) = / k{g,q)U{q)dq. 

qeSE(3)SE(3) SE{3) 

1. Operator <P is legal, so from the first identity in Definition [1] we have 

Vg'eSB(3)V[>gL«(SiS(3)) : (.^oCgO u){g') = {Cg o 3 o U){g'). (9) 

This holds for all U and by writing m in integral form, one can check that 
k{g~^g',q) = k{g',gq), k{e,g-^q) = k{g,q), ioi all g,g',q G SE{3). (10) 

2. The second identity in Definition [U together with the unimodularity of the 
group SE{3) gives k{g,qh~^) = k{g,q) = k{gh',q), so k must be right- 
invariant under subgroup El with respect to both entries. 

3. Suppose U G^tj, then U{y',Rn’) = U{y',n') and 


(#(17))(y,R„)= / K(y,R„,y',R')C/(y',R')dy'dR' 

SEi3) 

= 2tt f K(y,R„,y',R„/)17(y',i?„/)dy'dcr(n'). 

R3xiS2 


( 11 ) 


Now p(g) = K{g, e), p{g) = K{e,g). Left-invariance then implies ([5]). □ 

Corollary 1. The exact scale space kernels pt ■ SE{3) -G andpt : —>■ 

K.'*" satisfy the following symmetries 


Pt(y, n) =pt(Re,,ay, Re,.an) andpt(/i ^g) = pt{g) = pt{gh') (12) 
for all t > 0, a > 0, (y, n) S xi S'^, g G SE(3), h, h' G H . Moreover, 

Pt{y, n) = pt{-Rly, R^e^) and pt{g) = Pt{g~^)- (13) 


Proof. The symmetry (fT^ is due to Lemma [TJ The second symmetry (IT!?!) is due 
to refiectional invariance (^ 13 ,^ 14 ,^ 15 ) >->■ (—yla, —yl 4 , —yls) in the diffusion (|1]) 
and reflection on the Lie algebra corresponds to inversion on the group. □ 

Remark 1. Note that (IT!?!) in terms of k would be: fc(y, n, 0, e^) = fc(0, e^, y, n), 
by the relation k{y, n, y', n') = p(R^,(y — y'), R^,n). This means that k defines 
a symmetric measure: evaluation in (y, n) of a kernel centered around the unity 
element (O.e^) should be equal to evaluation in the unity element of a kernel 
centered around (y, n). 
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Now the Gaussian kernel approximation, based on the logarithm on SE{3) and 
the theory of coercive weighted operators on Lie groups HZ], presented in earlier 
work |ini ch: 8 ,thm.ll] does not satisfy this symmetry. Next we will improve it 
by a new practical analytic approximation which does satisfy the property. 


2.3 Ne-w vs. Previous Kernel Approximations 

We first present two existing approximations for pt and pt , that do not automat¬ 
ically carry over the properties we have shown in the previous section. A possible 
approximation kernel for pt is based on a direct product of two 5'£’(2)-kernels, 
see m Eq.(lO), Eq.(ll)], to which we refer as This approximation is 

easy to use since it is defined in terms of the Euler angles of the corresponding 
orientations. However, in Eig. [3]we show that the symmetries described before 
are not preserved hy and errors tend to be larger when D 44 and t increase. 

Therefore we move to an approximation for the S'i?(3)-kernel pt, for which we 
show that it can be adapted such that it has the important symmetries. We 
need the logarithm on SE{5) for this approximation, so first consider an ex¬ 
ponential curve in SE{3), given by y{t) = 50 exp ^ EA^. The logarithm 
log 5 £;( 3 ) : SE{3) Tf.{SE{3)) is bijective, and it is given by 

exp c'{g)A^ ) ^ ^ (14) 

and we can relate to this the vector of coefficients c{g) = (g), (g)) = 

{c^{g), ■ • ■ ,c^{g))^. We define matrix f2 as follows: 

/ 0 -c® c5 \ 

c® 0 -c^ , R = e*^, R e 5'Cl(3). (15) 

\-c® c-* 0 / 


logS£;(3)(ff) — logS£;(3) 


We can write R in terms of Euler angles, R = Re^^-yRe ./jRe^.a- Let the matrix 
L? 7 ,/ 3 ,a be such that = R-ej, 7 R-e„./ 3 Re 2 .a- The spatial coefficients 

c(i) = ^ ^ are given by the following equation: 



where is the (Euclidean) norm of Then another approximation for 

the kernel pt is given by [5] : 




1 I *°gSE(3)(g)IC33 

033044)'^ ^ 


(17) 


with the smoothed variant of the weighted modulus, [TOl Eq.78,79], given by 


r.l|2 


r.2|2 


r.6|2 


D:i^D 


33-^44 


D, 


44 


(C3)2 


r.4|2 


D. 


-k 


-k c 


,512 


33 


E 


44 


■ (18) 


I logs£;(3)(-)|D33,D44 — 
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Now the difficulty lies in the fact that the function U is defined on the quotient for 
elements (y,Rn), where the choice for a in the rotation matrix Rn (mapping 
02 onto n € is not of importance. However, the logarithm is only well- 
defined on the group SE{3), not on the quotient xi and explicitly depends 
on the choice of a. It is therefore not straightforward to use this approximation 
kernel such that the invariance properties in Corollary [T] are preserved. In view of 
Corollary [U we need both left-invariance and right-invariance for pt{g) under the 
action of the subgroup H. As right-invariance is naturally included viapt(y, R) = 
Ptiy, Rgz), left-invariance is equivalent to inversion invariance. In previous work 
the choice of a = 0 is taken, giving rise to the approximation 

PP’'’^(y.n(/3,7)) =Pt°®(c(y,R.e,,7R-ej„/3))- (19) 

However, this section is not invariant under inversion, as is pointed out in Fig. [31 
In contrast, we propose to take the section a = — 7 , which is invariant under 
inversion (since (Re,j, 7 Re„,/ 3 Re 2 ,- 7 )~^ = R-e;,. 7 R-e;,,-;SR-ej,,- 7 )- Moreover, this 
choice for estimating the kernels in the group is natural, as it provides the weak¬ 
est upper bound kernel since by direct computation one has a = —7 =X ce = 0 . 
Finally, this choice indeed provides us the correct symmetry for the Gaussian 
approximation of ptiy, n) as stated in the following theorem. 

Theorem 1. When the approximate kernel on the quotient is related to 

the approximate kernel on the group p\°^ by 

pJ‘®“(y,n(/3,7)) := p'°®(c(y, Re,,7Re„./3Re,-7)), (20) 

i.e. we make the choice a = — 7 , we have the desired a-left-invariant property 
pr“(y,n)=pr“'(Re..a'y,Re.,a'n), a'e[0,2^]. (21) 


and the symmetry property 


pr“(y,n)=pr“'(-Rny,Rne.) 


( 22 ) 


Proof. We start by proving the a-invariance. Following definition (j20D we have 


^new 

Pt 


(^ 62 : ,a'Y 1 


■Gz 


n(/3,7)) =:pl°®(c(R.e 


1^62: , 7 +a' 




—( 7 +Q:' 


))) 


Recall that the matrix is defined by R.^^; 3 _c( := R-e 2 , 7 R-e„,/ 3 R-e 2 .a = 

Therefore in our case, where we choose a = — 7 , we find for all a': 


R-e„,/3 — Re;j, 7 +Q:'® 


i,-h+<.')R^ 


■e,,7+a' = e 

_ R-7,(a 


7 + q ;',/3, —(■7+ q :') 

7 + q :',/3, —(7 + 0:') 


)R. 


)R. 


z .T + Q 

Re.. 


We see that Re^ 


(17.y_/3^_.y)Rg - From this we deduce: 


( 2 ) _ T} (2) 

^ 7 +a',/ 3 , —( 7 +ck') 7 ’ 


and g^^Q^f ^ 7 ,/ 3,—7 ■ 


(23) 
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Fig. 3. Top: 2 positions and orientations, rotated around the 2 -axis. The graphs display 
Pt(R.e 2 .c«'yi R-e^.o'n), for jjj green, in orange and p"®™ in blue. Parameters 

D 33 = 1, D 44 = 0.02, t = 4. Bottom: k{y, n, 0 , e^) should be equal to fc( 0 , 62 , y, n), see 
Remark[T] The absolute sum of the difference of the two evaluations on a 5 x 5 x 5 grid 
with at each point 42 uniformly distributed orientations is shown for j^prev ,2 

and fc"®™, corresponding to respectively, with coloring as 

above. Parameter D 33 = 1 and t, D 44 are as shown. 


and together with (HH) it gives 7 -^ 0 ',/?,-( 7 -Ha') = Com¬ 

bining this with (l23l) gives 


^/x,7+a',/9, —( 7 +a'') 


= Zi 


-X,7,/3,—75 


with Za> 



(24) 


It follows immediately that cf -I-C2, C3, C4-I-CI and cq are independent of a'. The 
proof for a-invariance is completed by stating that given a' G [ 0 , 27r]: 

p7“'(Re,.,a'y, Re,.a'n(/3, 7 )) = (log 5 £;( 3 ) (R-e,. .a'N, R 7 +a',/ 3 ,-( 7 +a') )) = 

Pt°®(^J'logs_E(3)(y,R7,/3,-7)) =Pt°®(log5£;(3)(y,R7,/3,-7)) =pr”'(y,n(/3,7)). 

The symmetry property (1221) directly follows from the fact that 


9 


-1 


6XPSB(3) 



(y,R)-i = (-R-V,R^'), 


the fact that our weighted modulus on Te{SE{3)) is invariant under reflection 
c I—H (—c), and the fact that the section a = —7 is invariant under inversion. □ 
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Fig. 4. Blue frame: Glyph field visualization, recall Fig. [T] of the FOD obtained with 
CSD from an artificial dataset (top left) can be improved significantly with enhance¬ 
ments (top right). Especially alignment on the boundary of the bundle is improved. This 
results in a more plausible tractography result (bottom). Here, D 33 = 1, D 44 = 0.02, 
t = 4. Red frame: example of the coherence quantification of a fiber bundle, the Optic 
Radiation. 


3 Neuroimaging Applications 

The newly proposed enhancement kernel is used within the pipeline depicted in 
Fig-Efor further processing of the dMRI data. There are two places where this is 
useful: for enhancement of the FOD obtained with e.g. CSD, and for quantifying 
the coherence of a fiber within a bundle. 

a. Enhancing FODs: To illustrate the use of enhancement of FODs, we 
use the artifical IEEE ISBI 2013 Reconstruction Challenge dataset [TH]. CSD is 
applied to a simulated dMRI signal, yielding an FOD as in Fig. The trac¬ 
tography was done with MRtrix [19j . with seed points randomly chosen in the 
bundle. Only the fibers that pass through the indicated spheres are shown. The 
alignment of glyphs improves due to the enhancements and this results in a bet¬ 
ter tracking: bundles are reconstructed fuller and less fibers take a wrong exit 
from the bundles. 

b. Coherence quantification of tractographies: The result of a prob¬ 
abilistic tractography is typically less well-structured than in the deterministic 
case, recall Fig.E We employ our new approximation kernel p"®™ to construct an 

X scale space representation of the set of fibers from a tractography. This 
density is locally low for a fiber that (partly) deviates in position or orientation 
from the other fibers in the bundle, as we explain next. 
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Suppose we have a set of N streamlines (fibers), the fiber having Ni points 
yj € and we set the orientation nj = — yl)/\\ylj^i — Yi || G 5^. Next 

we use the set of streamline points |(y|, ±n^) | j = 1 ,... f = 1,..., as 
follows in the initial condition of our scale space evolution ([5]): 


U{y,n) 


1 


2 N Ni 


<7—1 i—1 j — 1 


(25) 


with Ntot the total number of fiber points and d(y_n) a d-distribution in x 
centered around (y, n). The solution is given by the convolution W{-,t) = 
(p"®*" *R 3 xi 52 U){-). Intuitively, letting this sum of ^-distributions diffuse ac¬ 
cording to (|ni) results in a density indicating for each position and orientation 
how well it is aligned, in the coupled R^ x ^^-sense, with the fiber bundle. 

In practice, we evaluate the convolution (p"®“ =i=r 3,<|52 U) only in the fiber 
points, giving kernel evaluations. However, this is reduced to evalu¬ 

ations, since the order of arguments in fc"®*" is irrelevant for every pair of fiber 
points, thanks to the included symmetry in the kernel. This implies a reduction 
of the computation time. Summation over one fiber then gives a value, which is 
the final measure for the fiber to bundle coherence (FBC). Fibers with low FBC 
can be removed from the tractography result and coherent fibers remain. A local 
FBC can be computed by summing over parts of the fibers instead of the entire 
fiber. Fig. |4] shows this local FBC for a tractography of the Optic Radiation, a 
brain white matter bundle connecting the Lateral Geniculate Nucleus and the 
visual cortex. It can be seen that where fibers deviate from the bundle, the FBC 
is lower and hence they could be excluded. 


4 Conclusion 


We have introduced a new approximation for the kernel of a linear scale space 
PDF (|ni) on R^ X S'^. This is done by the embedding of R^ x into SE{3), 
on which the PDF coincides with the Fokker-Planck equation for a Brownian 
motion process. Due to this embedding, the kernel is subject to two constraints, 
recall Corollary [TJ We have shown the application of the kernel in two dMRI 
neuroimaging applications. First, we have used the kernel implementation in 
combination with CSD to enhance fiber structures in dMRI and thereby improve 
tractography results. Secondly, the kernel was used to construct a measure for 
coherence of fibers in a fiber bundle. Thanks to the appropriate symmetries of the 
scale space kernels, this measure is based on symmetric distances on R^ x and 
a reduction of computation time is obtained. In future work we aim to quantify 
the improvement in tractography due to the enhancements, and we pursue the 
use of the FBC for better construction of the structural connectome |20) . 
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